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Abstract. - Maximum-likelihood methods are applied to the problem of absorption tomog- 
raphy. The reconstruction is done with the help of an iterative algorithm. We show how the 
statistics of the illuminating beam can be incorporated into the reconstruction. The proposed 
reconstruction method can be considered as a useful alternative in the extreme cases where the 
standard ill-posed direct-inversion methods fail. 



Introduction. The standard reconstruction method in present computerized tomo- 
graphic (CT) imaging is the filtered back-projection (FBP) algorithm which is based on the 
Radon transformation [[l]. Unfortunately FBP fails in case of missing projections and/or if 
strong statistical fluctuations of the counting numbers are present in the small detector pix- 
els. The latter situation occurs e.g. in neutron tomography if monochromatic neutron 
beams are applied in order to avoid beam artifacts || or at the investigation of strong absorb- 
ing materials. The cases of missing projections and incomplete data sets for monochromatic 
neutron beams have been already investigated in the past in detail by means of algebraic re- 
construction technique . Scattering data from a double crystal diffractometer have been 
used to reconstruct 2D scattering pattern and the results were compared with the standard 
FBP. With this algebraic approach one could reconstruct 2D pattern in spite of the lack of 
nearly 90 degrees of the scanning angle, whereas in such cases the FBP method entirely failed. 
The computing time, however, was extremely long (up to several hours), so that this method 
is useful for rather small 2D arrays (100 x 100 pixels) only. 

The new reconstruction method proposed in this paper can improve several tomographic 
applications in neutron optics which in many cases are limited by the weak intensity and the 
poor detector resolution. The use of well collimated pencil beams which are scanned across 
the sample surface could dramatically enhance the spatial image resolution but this method 
is only rarely used due the long measurement times Jt0[ . An improved reconstruction method 
can encourage new applications in neutron optics which often suffer from the low counting 
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Fig. 1 - Geometry of the experimental setup 
Fig. 2 - Definition of coefficients dj 



numbers. Generally the new algorithm can achieve better reconstruction results or reduce the 
scanning time in neutron optics and in medical and biological CT imaging. 

LinPos tomography. Basic notions and the geometry of experimental setup are as 
follows. Let us assume that the sample is illuminated by parallel monochromatic pencil 
beams, see Fig. [lj Data consist of the number of particles counted behind the sample for M 
different scans - each scan being characterized by horizontal position h and rotation angle ip. 
Alternatively, a broad illuminating beam combined with a position-sensitive detector (CCD 
camera) placed behind the sample can be used. In that case h labels pixels of the camera. For 
the sake of simplicity a collective index j = {h, ip} will be used, hereafter, to label the scans. 

Mean number fij of particles (intensity) registered in j-th scan is given by the exponential 

law 

hj = n exp(- J n(x,y)dsj), (1) 

where fio is the intensity of the incoming beam, /j,(x, y) is the absorption index (cross section) 
of the sample in position {x, y}, and the integration is the path integration along the pencil 
beam. This exponential attenuation law is a good approximation if scattering can be neglected. 
The beam hardening artifacts would also modify Eq. ([!]) but this complication can be avoided 
experimentally by the use of monochromatic beams ||. For practical purposes, it is convenient 
to discretize Eq. (fy) as follows, 



JV 

lj = n exp(- ^ MiCij; 

i=0 



(2) 



The sample is now represented by a 2D mesh. Each cell is assumed to have a constant 
absorption index. The variables are now N numbers fa specifying absorption indices of those 
cells. Overlaps between beams and cells are stored in the array {cy}, see Fig. pi 

Let us first ignore the statistics of the illuminating beam, and assume that the counted 



numbers of particles {rij} do not fluctuate, 



Vj. Taking logarithms of both sides 



of Eq. (|2j), one obtains a system of M linear algebraic equations for N unknown absorption 
coefficients fa: 

fj=Pj, i = l...M, (3) 
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where we defined, 



fj = - m — , pj = ■ ( 4 ) 



Notice that problem (Q) is a linear and positive (LinPos) problem. Positivity follows from 
the fact that no new particles are created in the sample. Although direct inversion of Eq. (|J) 
is possible for N > M, the solution is not always positively defined. A negative value of a 
reconstructed Hi would suggest that particles were being created in the i-th cell in the course 
of the experiment, which would obviously be a wrong conjecture. This problem can be avoided 
if the problem (|J) is solved in the sense of maximum likelihood (ML) on the space of physically 
allowed absorption coefficients. In this approach one considers the data / and the prediction 
of the theory p as two probability distributions. One looks for absorption coefficients {/!;} 
that minimize the Kullback-Leibler "distance" 

d(/,p) = -E^ ln T ( 5 ) 

j h 

between the data / and the theory p. Here a little extra care is needed since p and / are 
generally not normalized to unity. The minimum of the Kullback-Leibler distance corresponds 
to the maximum of the maximum likelihood (ML) functional |U| 



p 3 " fi 



J2kPk 



(6) 



that quantifies the likelihood of the given distribution {/^} in view of the registered data. We 
seek the maximum-likely distribution of the absorption indices. A convenient way how to find 
it is the celebrated Expectation-Maximization (EM) iterative algorithm |l2[^3|. 

= H(/i w ) • n n , (7) 

where 

i' 

and fx is some initial strictly positive distribution (if 1 ' > 0, i = 1 . . . N. A nice feature of EM 
algorithm is that its convergence is guaranteed for any input data fj |l1|| . For this reason it 
became a valuable tool in many inverse problems which can be reduced to the form of Eq. (^|), 
e.g. in positron emission tomography |p^-[l6|]. The original derivation of EM algorithm is 
based on alternating projections on specially chosen convex sets of vectors. However, one 
could directly use the calculus of variations to derive the necessary condition for the extreme 
of the functional (^). Iterating these, one eventually arrives at the EM algorithm again. An 
advantage of this alternative derivation is that it can be also applied to more realistic physical 
models of the actual absorption experiment. One such possible generalization will be studied 
in the following section. 

Tomography with Poissonian signals. - Real signals are not composed of a sharp number 
of particles. For instance, two signals often used in experiments — beam of thermal neu- 
trons and laser light — both exhibit Poissonian fluctuations in the number of particles. Also 
monochromatic neutron beams are correctly described by Poissonian statistics if the detected 
count events occur mutually independently fll7| |. The knowledge of the true character of sig- 
nal illuminating the sample is a useful piece of prior information, which can be utilized for 
improving the performance of ML tomography. 
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As the Poissonian character of the signal is preserved by the process of attenuation, the 
counted numbers of particles behind the sample are random Poissonian variables. The corre- 
sponding likelihood functional reads, 



nfe-**. (9) 

This is the joint probability of counting {rij} particles. Mean values {rij} obey the exponential 
law (0) as before. They depend on the absorption in the sample {fij} that is to be inferred 
from the data. The necessary condition for the extreme of the likelihood (^) can be derived 
using the calculus of variations. The extremal equation can be shown to have the same vector 
form as the extremal equation of the LinPos problem (Q) . The vector R now becomes 

^(Poisson) = riQ ^2 c .. exp (_ ^ IH'Ci'j). (10) 

j' 3 

When the input intensity no is not known, it can be estimated together with the absorption 
of the sample: 

n o = ^ f—. r- (11) 

j i 

Poissonian tomography is intrinsically a nonlinear problem. This has serious consequences for 
the convergence properties of the iterative algorithm (Jfy and (fTcf). Instead of converging to 
a stationary point it might end up in oscillations. Typically such convergence problems arise 
in the presence of very noisy data. When this happens one can always decrease the length of 
the iteration step as follows: Ri — > Rf, i = 1 . . . M, < a < 1. Of course, any solution to 
the regularized problem is also a solution to the original problem. 

Discussion. - Generally, the reconstructed image will depend on which ML method is 
chosen to process the data; see the apparent difference between Eqs. (||) and (|l0|). It is inter- 
esting to look more closely at the origin of this difference. Consider a tomographic setup with 
a Poissonian beam. Then the Poissonian algorithm should provide a better reconstruction 
than the LinPos algorithm which have been derived under the assumption of non-fluctuating 
signals. The LinPos reconstruction consists in minimizing the Kullback-Leibler distance be- 
tween the data / and theory p. When logarithms of the counted numbers of particles are 
chosen to be the input data rather than counted data itself, one arrives at the EM algorithm 
(0) and (||). Taking logarithms of actual data makes the problem linear and considerably 
simplifies the reconstruction. However, one could, instead, directly minimize the Kullback- 
Leibler distance between the counted data rij and the corresponding theory p'j = no exp(-pj). 
Interestingly enough, the extremal equations associated with this variational problem are the 
same as Eqs. (Q) and (|l(]) derived above from the Poissonian theory @. Choosing rij instead 
of fj as the data is equivalent to taking the Poissonian statistics of the signal into account! 
The difference between the LinPos and Poissonian ML reconstructions can thus be traced 
down to whether the measured data are used directly or not. Tampering with data prior 
to reconstruction may speed up and facilitate the whole process of reconstruction but some 
information about the object might get lost. 

Comparison with standard methods. - In a real experiment there are many factors that 
could influence the quality of the measured data and therefore also on the result of the tomog- 
raphy. Misalignments present in the experimental setup, instability of the illuminating beam, 
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Fig. 3 - The object. 

white spots and damaged detector pixels can be such factors, to name a few. To avoid this 
problem we replaced the experiment by a simulation. The data were generated on a computer. 
The artificial object used in the simulation is shown in Fig. [| The object is a circle made of 
a homogeneous material with many small round holes drilled through it. One additional rect- 
angular piece of material was removed from the circle to make it less symmetric. Absorption 
index of the material was chosen in such a way that the maximum attenuation along a beam 
was close to 50% of the input intensity. 

In the simulation, the object was subject to five different experiments. Their parameters 
are summarized in Table || First four experiments correspond to the ideal situation of a 
very high beam intensity where the Poissonian detection noise can safely be ignored. The 
last reconstruction simulates more realistic conditions with 2000 counts per pixel in the open 
beam. Notice that a relatively small number of rotations is chosen for all five experiments. In 
this regime the Radon transformation is expected to yield bad results and the improvement of 
the maximum-likelihood tomography upon the standard technique should be most prominent. 
This regime is also important from the practical point of view. Doing more rotations implies 
a longer measurement time and more radiation absorbed by a sample. The latter may be 
an important factor if the imaging of biological samples is considered. So, imaging costs and 
damage done to a sample due to radiation might be reduced provided the improvement of the 
reconstruction technique gives comparable resolution with less data. 

Reconstructions from the simulated data are shown in Figs. || and ||. The simulated 
data were first processed using the IDL imaging software (Research Systems Inc.) which 
implements the standard FBP algorithm (Radon transform), see Fig. ||. This software is 
one of the industrial standards in the computer assisted tomography. The same data were 
then processed using our iterative algorithm based on the maximization of the Poissonian 
likelihood function, see Fig. 0. In the absence of noise, see cases (a)-(d), the fidelity of a 
reconstruction depends on two main factors — the spatial resolution of the detector, and the 
number of rotations used. It is apparent from Figs. |] and |^ that the latter factor is more 

reconstruction angles pixels intensity 



a 


13 


161 


oo 


b 


19 


101 


oo 


c 


20 


101 


oo 


d 


7 


301 


oo 


c 


15 


161 


2000 



Table I - Quality of the input data. The last column shows the mean number of counted particles 
per pixel in the incident beam. 



6 



EUROPHYSICS LETTERS 




important of the two. Very small number of angles cannot be compensated by an increased 
spatial resolution of the detector, compare e.g. cases (c) and (d), and reconstruction (d) is by 
far the worst one. However, ML tomography is much less sensitive to the number of angles 
than the standard filtered back-projection. Even the large rectangular hole in the object 
is hardly perceptible in Fig. whereas it nicely shows in the ML reconstruction Fig. ||d. 
ML reconstructions are superior to the standard ones also in cases (a)-(c); notice that the 
reconstruction Fig. ||c done with as few as 20 different angles is nearly perfect. 

Benefits of the ML tomography are fully revealed when the detected data are noisy. This 
is case (e) in Tab. |. Standard filtered back-projection applied to noisy data faces serious 
difficulties. This is due to ill-poseness of the Radon transformation where data are integrated 
with a singular filter function. Obviously such deconvolution greatly amplifies any noise 
present in the data. Having little or no prior information about the object it is difficult 
to tell true details of the object from artifacts. ML tomography gives much better results. 
Since noises are incorporated into the algorithm in a natural and statistically correct way 
artificial smoothing is not needed. Notice in Fig. ||e that noisy data yield a little distorted 
but otherwise clear image unlike the corresponding very noisy standard reconstruction shown 
m Fig. §3. This is a nice feature of the intrinsically nonlinear ML algorithm which, in the 
course of reconstruction, self-adapts to the registered data and always selects the most likely 
configuration. 

Finally let us emphasize that apart from the size of the reconstruction mesh N |Tsf| there 
are no free parameters left in the ML algorithm to play with. This prevents one from in- 
terfering when the reconstructed image "looks bad." This also makes the whole procedure 
more objective, which is a necessary presumption for the investigation of ultimate limits of 
reconstruction schemes. 

Conclusion. We presented a new reconstruction method for CT imaging based on 
the iterative maximization of the Poissonian likelihood. For small number of scans and/or 
short measurement time this method was shown to yield a significant improvement upon the 
standard filtered back-projection algorithm. This could be important for CT imaging with low- 
intensity beams, and for applications where strong irradiation of a sample during the scanning 
should be avoided. One area where reconstruction techniques of the type discussed in this 
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paper would be very useful are coherent reconstruction techniques such as interferometric 
phase tomography with X-rays [|l9| or neutrons [^lj , or neutron holography |^] . There is 
hopefully more to come. 
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